Development of a two-particle self-consistent method for multi-orbital systems and its 

application to unconventional superconductors 



CO 

o 

CN 
C 

in 



Hideyuki Miyahara 1 , Ryotaro Arita 1,2 , and Hiroaki Ikeda 3 
1 Department of Applied Physics, University of Tokyo, 
7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan 
2 JST PRESTO, Kawaguchi, Saitama, 332-0012, Japan and 
3 Department of Physics, Kyoto University, Kyoto 606-8502, Japan 
(Dated: January 16, 2013) 

We extend the two-particle self-consistent method proposed by Vilk and Tremblay (J. Phys. I 
France 7, 1309-1368 (1997)) to study superconductivity in multi-orbital systems. Starting with the 
sum rules for the spin and charge susceptibilities, we derive self-consistent equations to determine the 
renormalized effective interactions. We apply this method to the two-orbital d x i _ y n-d- iz n _ r i model 
for La2Cu04 and the five-orbital d-model for LaFeAsO. Comparing the results with those of the 
random phase approximation or the fluctuation exchange approximation in which vertex corrections 
are ignored, we discuss how the vertex corrections affect the pairing instability of La2Cu04 and the 
dominant pairing symmetry of LaFeAsO. 
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I. INTRODUCTION 

Since the seminal studies by Suhl 1 and Kondo, 2 super- 
conductivity in multi-orbital systems has been one of the 
major topics in condensed matter physics. So far, many 
kinds of multi-orbital superconductors such as MgB 2 
(Ref. 3), alkali-doped C 60 (Ref. 4), Na x Co0 2 • 2/H 2 
(Ref. 5), Sr 2 Ru04 (Ref. 6), iron-based superconductors, 7 
and heavy fermion superconductors 8 have been discov- 
ered and studied extensively. Theoretically, a variety 
of exotic unconventional pairing mechanisms going be- 
yond the Migdal-Eliashberg theory 9 have been proposed. 
For example, it has been considered for the cobaltate su- 
perconductor that the Hund's coupling (which of course 
does not exist for single-orbital systems) induces triplet 
superconductivity, 10 and it has recently become an is- 
sue of hot debates whether orbital fluctuations mediate 
superconductivity in the iron-based superconductors. 11 

To investigate these fascinating possibilities, accurate 
calculations of superconductivity in correlated multi- 
orbital models are indispensable. Among many available 
approaches, from the weak coupling side, one often starts 
with the random phase approximation (RPA). Since the 
pioneering work for the single-band Hubbard model by 
Scalapino et al., 12 RPA has been successfully applied to 
various multi-orbital systems. The fluctuation exchange 
approximation (FLEX) developed by Bickers et al., 13 
which includes the self-energy correction self-consistcntly, 
has been also widely used. Here, the self-energy is calcu- 
lated in the manner of Baym and Kadanoff, 14 and conser- 
vation laws for one-particle quantities such as the total 
energy and momentum are satisfied. However, due to 
the absence of vertex corrections, FLEX violates conser- 
vation laws for two-particle quantities. 

Recently, several diagrammatic methods, which take 
into account some vertex corrections, have been 
proposed. 15 ' 16 Among them, the two-particle self- 
consistent method (TPSC) proposed by Vilk and 
Tremblay 17 is a promising approach in that it is com- 
patible with conservation laws in the two-particle level. 



In this method, vertex corrections in the charge and 
spin channel are assumed to be momentum and fre- 
quency independent, and they are determined in such a 
way that the correlation functions meet their sum rules. 
With this numerically inexpensive treatment, it has been 
demonstrated for the single-band Hubbard model that 
TPSC shows good agreement with quantum Monte Carlo 
(QMC) calculations. 

In this paper, we formulate TPSC for the multi-orbital 
Hubbard model. First, we derive a series of equa- 
tions to determine the vertex corrections in the spin 
and charge channel, and then apply this method to a 
two-orbital model for La 2 _ x (Sr/Ba) x Cu04 and a five- 
orbital model for F-doped LaFeAsO. Recently, the two- 
orbital model (which we call the d x 2_ y 2-d^ z 2_ r 2 model) 
was studied by FLEX 18 to give an insight into the ma- 
terial dependence of superconducting transition temper- 
ature (T c ). While FLEX successfully describes the dif- 
ference between La 2 _ r (Sr/Ba) a ;Cu04 (T c ~ 40 K) and 
HgBa 2 Cu04+a (T c ~ 90 K), it underestimates the pair- 
ing instability for La 2 _. c (Sr/Ba) a; Cu04 and T c is much 
lower than the experimental value. We show that, in the 
present multi-orbital TPSC calculation, the inter-orbital 
scattering enhances the d-wave instability and reason- 
able value of T c is obtained for the intermediate coupling 
regime. For the five-orbital model, it has been extensively 
studied by RPA 19 ' 20 ' 22 and FLEX. 21 There, strong spin 
fluctuation has been shown to mediate the s-wave su- 
perconductivity with sign changes (the so-called s±-wave 
pairing) . On the other hand, recently, it has been pointed 
out that vertex corrections can enhance orbital fluctua- 
tions, which mediate s-wave superconductivity without 
sign changes (the s++-wave pairing). 16 In this paper, we 
show that orbital fluctuations arc enhanced in TPSC, 
while the dominant pairing symmetry is still s± when 
the system resides in the weak coupling regime. 

This paper is organized as follows. In Sec. II, we for- 
mulate multi-orbital TPSC for the Hubbard model. We 
discuss how we calculate the charge (orbital) and spin 
correlation functions. In Sec. Ill, we show the results for 
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the two- and five-orbital Hubbard model and the sum- 
mary of the present study is given in Sec. IV. 



II. METHOD 
A. Model 

The Hamiltonian of the multi-orbital Hubbard model 
is given by 



H = ^e M n /1(7 (r)+ ^ C' c ^( r ) c ^( r ') 

+J2[ U Y1 n ^ ( r ) n *4 ( r ) + ^'EE n ^ ( r ) 

- J ]T M') ' S„(r) + J' 2 4 t ( r ) c L( r ) c ^( r ) c -t( r ) 



where cj^r) is a creation operator of an electron 
with spin a and orbital fi at site r, and n^ a {r) = 

4 CT (r) C/1(T (r), S M (r) = (cj^r), c^(r)) CT (c„ t (r), c^(r)) T 
with the Pauli matrices <x. The on-site Coulomb inter- 
actions, U, U', J, J' denote the intra-orbital, inter-orbital 
Coulomb repulsions, the Hund's exchange, and the pair- 
hopping term, respectively. 



B. Two-particle self-consistent method for the 
single-orbital Hubbard model 

Let us start with a review of TPSC for the 
single-orbital Hubbard model formulated by Vilk and 
Tremblay. 17 The central quantities in this method are 
the spin and charge correlation functions. In the non- 
magnetic state, the system holds SU(2) symmetry, and 
the spin-spin correlation functions do not depend on the 
spin directions. Thus we consider the z component of 
the spin operator S z (r) = n-j-(r) — ni(r) and the charge 
operator n(r) — nf(r) + n^(r). In RPA, the spin and 
charge correlation functions are evaluated as follows: 



2X°(<?) 



2X°(9) 



Xrpa(9)- i^ Ux o {q y Xrpa(s)- i + U X °(qy W 
with the irreducible susceptibility, 



T 



(«) = t Eg°Wg°(h ? ), 



where T and N are temperature and number of sites in 
the system, and G°(k) = l/(ie n + fi — e(k)) is the bare 
Green's function with chemical potential /i and energy 
dispersion e(k). Here, we have introduced the abbrevi- 
ations k = (k, ie n ) and q — (q, iu n ) being the fermionic 
and bosonic Matsubara frequencies, respectively. 



Note that the RPA violates the Pauli principles, and 
it does not fulfill the following two sum rules, 



T 
N 



T 



£x sp (9) = (Mr) -%(r))(n t (r) - n ; (r))) 
i 

= n-2(n t n 4 ), (2a) 



X £ ^ Ch (?) = «"t« + %(r))(n t (r) + H (v))) n 2 



n + 2(n t n i ) - 



(2b) 



which are exact relations derived via the Pauli princi- 
ples, (n a (r) 2 ) = (n a (r)) (see Appendix A). Here, n is 
the particle number per site and for nonmagnetic states, 
(n,f(r)) = (nj.(r)) = n/2. Note that the double occu- 
pancy, (rif(r)n_i(r)) = (n-j-nj,) is also translation invariant 
and does not depend on site r. 

In TPSC, to meet the above conditions [Eqs. (2)], we 
introduce two independent effective interactions, U sp for 
the spin channel and U ch for the charge channel. Then 
the full susceptibilities of Eq. (1) are replaced with 



X sp (q) 



2x°(g) 



X ch (q) = 



2 X °(q) 



l-U s Px°(q)' l + U ch X °(q)' 

Finally, we put the following ansatz: 



U sp 



u, 



(3) 



(4) 



which is compatible with the equations of motion (see 
Appendix B). Equations (2), (3), and (4) provide a set 
of self-consistent equations in TPSC. Namely, U ch , U sp , 
and (n-j-nj,) are self-consistently determined for given n 
and U. The one-particle Green's function and self-energy 
are calculated by 



G(k) = G°(k) + G°(k)E(k)G(k), 
1 T 
47V 



E W = E [ USP X SP W + U ch X ch (q)u] G(k - q) 



For the single-band Hubbard model, it has been demon- 
strated that TPSC agrees well with QMC. 17 



C. Extension to multi-orbital systems 

Let us here formulate TPSC for the multi-orbital Hub- 
bard model. Hereafter, we follow the matrix form em- 
ployed in Rcfs. 10 and 23. The irreducible susceptibility 
is defined as 



(5) 



which can be considered as a matrix element with a row 
X/2 and a column of a matrix X°(q)- I n the non- 
magnetic state, the system is invariant for spin rota- 
tion, and then 2x spz (q) = X sp± (q) holds, where % sp± (<z) 
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is the in-plane correlation function between (r) — 

SL(r) ± *S^(r) = 4 CT (r)c M5 (r) with a =| or |. In 
TPSC, similar to the single-orbital case, the spin and 
charge susceptibilities are given by 



X sp (q) = (1 - X°(g)U ,p )- 1 2 X °( 9 ), (6a) 

x ch (9) = (i + x°(9)u ch )- 1 2 X (g), (6b) 

where U sp ^ ch - ) is the renormalized effective interaction 
matrix for the spin (charge) channel. 10 ' 23 For two-orbital 
systems, for instance, these are represented as 



ccptibilitics and the double occupancy: 



T 



T 
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+ 2(n fi - r n fi i) - n 



(9a) 



E^^H^T ^ + c^c I/4 .)(4 t c A(t + c^c M )) 

9 

=<n Mt (l - ?m)) + ('VlX 1 - 

= | (E + W-^)) • (9b) 
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where Uf^ 



Q^uuuu) i g *be intra-orbital Coulomb in- 
teraction; Ufv„ v (Ufyuv) with /j / i/, the inter-orbital 
Coulomb interaction; J sp ( J ch ), the Hund's coupling. 10 ' 23 
In the present study, for simplicity, we ignore the Hund's 
coupling in the charge channel, namely, J ch = 0. 24 In 
RPA, one employs the unperturbed bare vertex as fol- 
lows, V$$ = U, U s $y = U^P = U', and J s ^ = 
J. 



Next let us consider the sum rule for multi-orbital sys- 
tems (see Appendix A). In the n-orbital Hubbard model, 
there are n 4 sum rules for x sp (<?) an( l X ch (<z)- Among 
them, we use the following equations to determine U sp : 



Finally, as in the single-band case, we introduce the 
following ansatz between the two-particle quantities and 
the interaction parameters (see Appendix C); 



j rsp _ (jWWj 

U%.....= ,^ n ^\ U', 



U sp - V 



{n<jn){n Sv ) 



(U'-J). 



(10a) 
(10b) 
(10c) 



Equations. (6)-(10) are a set of self-consistent equations 
in the multi-orbital case. 



D. Eliashberg equation 

Superconductivity has been studied by the following 
linearized Eliashberg equation, 



\A U r(k)= ^2 Vi mimi v{k,k')G mim2 {k') 



(11) 



T 



N 

T 

N 



N 



Ex^)= 2 <it c ^Mt) 

=2(n Mt ) - 2(n Mt n H ), 



(8b) 



E*w^('?H 2 ( ? Vt r M) - 2(n Alt n Ali ). (8c) 



Note that the intra-orbital component of the sum rule 
has the same form as that of the single-orbital Hubbard 
model. For the inter-orbital components, we use x sp± (<?) 
rather than X spz (l)j since they can be expressed in terms 
of the density operators. 

For U ch , we use the following sum rules for the charge 
susceptibilities which can be represented by the spin sus- 



Eigenstate A/;/ (fc) with the largest eigenvalue A was nu- 
merically evaluated by the power method. The supercon- 
ducting transition occurs at the temperature for which A 
becomes unity. Here, Gw (k) is the dressed Green's func- 
tion, 

G lv (k) = G%(k) + Gf m (k)E mm ,(k)G mrv (k), (12) 
and the self-energy (k) is given by 



1 T 

In 



(13) 



ch 



U c V h (<7)U c 



ImV 



Gmm'(k - q). 



In the present study, we omit the Hartree-Fock term, 
since a part of its contribution is already considered in 
the one-body part of the Hamiltonian, which is derived 
from density functional calculation. 
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The effective interaction Vu> mm >(k, k') for the spin- 
singlet pairing can be expressed in a matrix form as fol- 
lows: 



La2Cu04 



V(fc) = — u sp x sp (fc) u* p + -U ch x ch (fc)UJi h 



2 o 2 o > 

ych 



(14) 



where Uq P and Uq are the bare vertex in the spin and 
charge channel, respectively. 10,23 



III. RESULTS 

Let us move on to the application of the multi-orbital 
TPSC method to the effective models for La2CuC>4 and 
LaFeAsO. Using the technique of the maximally local- 
ized Wannier functions, 25 these models are derived from 
first-principles calculations. In the density-functional 
calculations, we employed the exchange correlation func- 
tional proposed by Perdew et al., 2e and the augmented 
plane wave and local orbital (APW+lo) method as imple- 
mented in the WIEN2K program. 27 We then constructed 
the Wannier functions for the d bands around the Fermi 
level, using the WIEN2Wannier (Ref. 28) and the wan- 
nier90 (Ref. 29) codes. 



A. La2CuC>4 

Recently, the two-orbital d x 2_ y 2-d 3z 2_ r 2 Hubbard 
model for the cuprates were studied to understand the 
material dependence of T c by FLEX. 18 There, the energy 
difference between the d x 2_ y 2 orbital and the d 3z 2_ r 2 or- 
bital was found to be a key parameter to characterize 
La2CuC>4 and HgBa2Cu04. Namely, in the FLEX calcu- 
lation for the two-orbital model, the pairing instability is 
stronger in the latter. On the other hand, in the former, 
the eigenvalue of the Eliashberg equation within FLEX 
does not reach unity down to T ~ 40 K. The purpose of 
this subsection is to examine how the vertex corrections 
in TPSC affect the superconductivity in La2Cu04. 

The band structure of the effective two-orbital model 
for La 2 Cu04 is shown in Fig. 1. We set U = 2.0 eV, 
U' = 1.6 eV, J = 0.2 eV, and n = 2.85. We employ 
64 x 64 fc-point meshes and 2048 Matsubara frequencies. 
Hereafter, orbitals 1 and 2 denote the d x 2_ y 2 and d 3z 2_ r 2 
orbitals, respectively. 30 

In Fig. 2, we plot temperature dependence of A, the 
maximum eigenvalue of the Eliashberg equation. While 
A does not show appreciable temperature dependence 
in FLEX, A is drastically enhanced at low temperature 
< 0.02 eV in TPSC. The characteristic enhancement of A 
in TPSC is attributed to the low-temperature behaviors 
of the spin and charge susceptibility. In Fig. 3, we plot 

Xmi(q^ = 0), Xm2(q> w = °)> and -Xmi(<li w = °) 
at T = 0.020 eV, which indicate that the system has 
a strong incommensurate spin correlation in the d x 2_ y 2 
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FIG. 1: (Color online) Band structure of the two-orbital 
model for La2Cu04. The model consists of the d x 2_ y 2 or- 
bital and the d 3 ,2_ r 2 orbital. The Fermi level is set at eV. 
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FIG. 2: (Color online) Temperature dependence of the max- 
imum eigenvalue of the linearized Eliashberg equation ob- 
tained by TPSC (red solid line) and FLEX (blue dotted line). 
U and J are 2.0 eV and 0.2 eV, respectively, and n is set to 
be 2.85. 



orbital and strong commensurate inter-orbital fluctua- 
tions. We here stress that there is no large peak in the 
charge susceptibilities in the RPA and FLEX calcula- 
tions, so that these enhanced inter-orbital charge fluc- 
tuations are purely due to the effects of vertex correc- 
tions. It should be noted that Xi22i(Q) has a negative 
peak around Q = (7r,7r), which works as attractive force 
between the two orbitals for d-wave pairing just like anti- 
ferromagnetic spin fluctuations. [Note that the spin and 
charge sectors in Eq. (14) have opposite signs.] As we 
will see below, there is a close correlation between the 
characteristic enhancement of A and the charge fluctua- 
tions. 

Figure 4 depicts temperature dependence of the max- 
imum value in the spin susceptibility, the orbital suscep- 
tibility, and its inverse. We can see that the enhance- 
ment of the peaks in X1212 an d — X1221 dominate over 
that of X1111 f° r T < 0.02 eV. This behavior comes 
from the fact that while U sp is renormalized substan- 



tially (*7ff n - 1.16 eV at T = 0.020 eV), (U$) 
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FIG. 3: (Color online) (a) Xiin(q,^ = 0), (b) Xi2i2(q,w = 0), and (c) -Xi22i(q,^ = 0) at T = 0.020 eV, where orbitals 1 and 
2 denote the d x a_ y a and d 3z 2_ r 2 orbitals, respectively. 
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FIG. 4: (Color online) Temperature dependence of the maximum value in (a) Xiiii; (b) 1X1212(1221) li an d (c) |Xi2i2(i22i) 



becomes ~ 1.55 eV, which is even larger than the bare 
value 1.2 eV used in RPA and FLEX. In fact, similar en- 
hancement of charge channel is also observed in the case 
of the single-orbital model. 17 

Figure 5 shows the gap functions for the d x i_ y i orbital 
and the d% z 2_ r 2 orbital (An and A22) at T = 0.022 eV. 
They have the d-wave symmetry, which is mediated by 
the dominant spin fluctuation, denoted by the black ar- 
row. This is the conventional situation where the orbital 
fluctuations remains small. 
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FIG. 5: (Color online) Gap function for the (a) d x 2_ y 2 or- 
bital and the (b) d 3z 2_ r 2 orbital at T = 0.022 eV. The gap 
functions have d-wave symmetry and the black arrow denotes 
the pair scattering mediated by antiferromagnetic spin fluc- 
tuations. 

In Fig. 6, we plot the gap functions at lower tempera- 



ture T = 0.018 eV. We see that characteristic structure 
emerges at (71", 0) and (0,7r) in the d 3z i_ r 2 gap function, 
due to the inter-orbital effective interaction, dominantly 
mediated by the orbital fluctuation XmiX?)) which con- 
nects An and A22 [Eq. (11)]. Since Xi22i(<?) takes a 
large negative value around (it, 7r), it cooperates with an- 
tiferromagnetic spin fluctuation to enhance the d-wave 
pairing instability. 
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FIG. 6: (Color online) Plots similar to Fig. 5 for T = 0.018 eV. 
The dashed black arrow denotes the pair scattering mediated 
by orbital fluctuations. 

This situation changes in the stronger coupling regime 
(U ~ 2.5 eV), where the system goes away from a su- 
perconducting instability. This is because the dominant 
spin/orbital fluctuations make the quasi-particle damp- 
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FIG. 7: (Color online) Temperature dependence of the maximum value in (a) Xiim O 5 ) IXi2i2(i22i) l> an< ^ ( c ) the maximum 
eigenvalue of the linearized Eliashberg equation obtained by TPSC for U — 2.5 eV, U' = 2.0 eV, and J = 0.25 eV. 



ing around (0, it) and (%, 0) significant, as is observed in 
the previous FLEX calculation. 18 Temperature depen- 
dence of the maximum value in the spin, orbital suscep- 
tibility and the maximum eigenvalue of the Eliashberg 
equation A for U = 2.5 eV, U' = 2.0 eV and J = 0.25 eV 
are shown in Figs. 7(a), 7(b), and 7(c), respectively. The 
maximum value in Xim always dominates over that of 
X1212 an d — X12121 ana - ^ does not show any enhancement. 



B. Iron-based superconductor: LaFeAsO 

Let us now apply multi-orbital TPSC to the iron-based 
superconductor, LaFeAsO. The recent discovery of high 
T c superconductivity in F-doped LaFeAsO 32 has stimu- 
lated a renewed interest in multi-orbital superconductors. 
As for the pairing mechanism of the iron-based supercon- 
ductors, several scenarios have been proposed. Among 
them, the possibility of the sign-reversing s±-wave su- 
perconductivity mediated by spin fluctuations 22 ' 33 have 
been extensively studied. While the s±-wave solution 
has been obtained in the RPA 19 ' 20 ' 22 or FLEX 21 calcu- 
lations for the five-orbital d- model, recently, it has been 
proposed that vertex corrections can enhance orbital fluc- 
tuations, and the s+_|_ -pairing without sign reversing be- 
comes dominant. 11 In this subsection, we discuss how 
vertex corrections in TPSC affects superconductivity in 
the five-orbital d- model for LaFeAsO. 

The band structure of the d-model is shown in Fig. 8. 
The bare coupling constants are set to be U = 1.5 eV, 
U' = 1.2 eV, and J = 0.15 eV. 34 The number of electrons 
n is 6.1. We employ 64 x 64 fc-point meshes and 2048 
Matsubara frequencies. Hereafter orbitals 1, 2, 3, 4, and 
5 denote the d 3z 2_ r 2, d xz , d yz , d x 2_ y 2, and d xy orbitals, 
respectively. 

To see that TPSC can give enhanced orbital fluctua- 
tions, we plot X2222(q> w = 0) and x^ 24 (q, w = 0) at T = 
0.015 eV in Figs. 9(a) and 9(b), respectively. Clearly, 
both susceptibilities have peaks around (ir,0) and (0, 7r), 
and the peak in the orbital susceptibility is higher than 
that of the spin susceptibility. Temperature dependence 
of these peaks are shown in Fig. 10. We see that the 
peak of X2424 (li w = 0) is more drastically enhanced 



LaFeAsO 




(0,0) (ti,0) (n,it) (0,0) 



FIG. 8: (Color online) Band structure of the five-orbital d- 
model for LaFeAsO. The Fermi level is set at eV. 



than that of x 2 222(q> w = °) for T > °- 01 cV - whilc 
X2424( ( li tL ' = 0) nas a broad maximum peak around 
T = 0.01 eV, X2222(qi w = 0) grows monotonously as 
temperature lowers. 

The enhancement in the orbital susceptibility comes 
from the vertex correction in the charge susceptibility. 
To make this point clear, in Fig. 11, we plot temperature 
dependence of the maximum value in X2A2a(Qi w = 0) and 
X2222 (Q; w = 0) obtained by RPA, for which the bare 
coupling constants are set to be U = 1.2 cV, U' = 0.96 
cV, and J = 0.12 cV. We see that while X2222( ( l' a; = 0) 
diverges around T = 0.01 eV, X2424( c l' a; = 0) has no 
significant temperature dependence. 

In Fig. 12, we show temperature dependence of the 
maximum eigenvalue of the Eliashberg equation. We see 
that the system has a superconducting transition around 
T ~ 0.005 eV. 

The associated eigenfunctions of the Eliashberg equa- 
tion at T — 0.015 eV (the gap functions) are shown 
in Fig. 13 for the three bands crossing the Fermi level. 
We see that these gap functions have the s± symme- 
try, indicating that the spin fluctuation is the primary 
glue of superconductivity. However, there is a notable 
difference between TPSC and RPA results in the am- 
plitudes of the gap functions on the Fermi surface. As 
we can see in Fig. 14, the gap amplitude is larger for 
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FIG. 9: (Color online) (a) X2222 ("^ ^ = 0) an d (b) X2424(*l> w = 0) at T = 0.015 eV, where orbitals 2 and 4 denote the d xz and 
d x 2_ y 2 orbitals, respectively. 
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FIG. 10: (Color online) Temperature dependence of the max- 
imum value in X2222 ( re d solid line) and X2424 (blue dotted 
line) obtained by TPSC. 
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FIG. 11: (Color online) Temperature dependence of the max- 
imum value in X2222 ( re d solid line) and X2424 (blue dotted 
line) obtained by RPA. The bare coupling constants are set 
to be U = 1.2 eV, U' = 0.96 eV, and J = 0.12 eV. 



RPA than TPSC. This indicates that there is a frus- 
tration between the orbital-fluctuation-mediated pair- 
ing and the spin-fluctuation-mcdiatcd pairing. Indeed, 
in TPSC, if we drop the contribution of the charge 
channel in the pairing interaction, namely consider only 



FIG. 12: (Color online) Temperature dependence of the max- 
imum eigenvalue of the linearized Eliashberg equation. The 
blue dotted line is a guide to the eye. 



V(fc) = -§U s Px sp (fc) Up p in the Eliashberg equation, 
then we find that the gap amplitude becomes large on 
the Fermi surface as in the RPA result [see Fig. 14(d)]. 



IV. SUMMARY 

To summarize, we have developed the two-particle self 
consistent method (TPSC) for the multi-orbital Hubbard 
model. We derived self-consistent equations to determine 
vertex corrections in the spin and charge (orbital) suscep- 
tibilities. We applied this method to the effective mod- 
els for La2Cu04 and LaFeAsO. We solved the linearized 
Eliashberg equation and found that vertex corrections 
play a crucial role in the multi-orbital superconductors. 

In the two-orbital d x 2_ y 2-d 3z 2_ r 2 model for La2Cu04, 
while FLEX shows much lower T c than its experimental 
value ~ 40 K, the present TPSC can increase T c dra- 
matically due to enhanced orbital fluctuations via vertex 
corrections for intermediate U ~ 2.0 eV. 

In the iron-based superconductor LaFeAsO, we have 
studied whether orbital fluctuations can be enhanced and 




FIG. 13: (Color online) Gap functions obtained by TPSC for the bands with the (a) second, (b) third, and (c) fourth Kohn- 
Sham energy at T — 0.015 eV. The black line and dotted green line represent the Fermi surface and nodes of gap functions, 
respectively. 



(a) 



71 



-71 



J 


W FS4^ 


) 


© c 




FS2 v 







-71 





9 X 



71 



o 
< 

C 

o 
s 

eg 

a 




Angle 



FIG. 14: (Color online) Gap function on the (a) Fermi surfaces by (b) TPSC, (c) RPA, and (d) TPSC without charge fluctuation, 
where 6 is the rotation angle from the k y axis. 



induce the s ++ -wave pairing within TPSC. Indeed we 
have found that some kinds of orbital fluctuations are 
enhanced by considering vertex correction and become 
even stronger than spin fluctuations. However, their or- 
bital fluctuations are not strong enough to cause the s++- 
wave pairing, and the obtained gap function has the s± 
symmetry, although the gap magnitude is relatively sup- 
pressed due to a frustration between two kinds of pairing 
interactions mediated by spin and orbital fluctuations. 
It is an interesting problem in future research whether 
the pairing symmetry changes for larger interaction pa- 
rameters. Another important future issue is a systematic 
comparison between the present multi-orbital TPSC and 
other (diagrammatic) methods which consider the vertex 
corrections. Finally, we stress the importance of vertex 



corrections in multi-orbital systems for cooperative and 
competitive phenomena between spin and orbital degrees 
of freedoms. 
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Appendix A: Definition of correlation functions 

In the single-orbital case, correlation functions for spin 
5 2 (r) and charge n(r) are defined as 

X sp (l,2) = (T T ^(l)^(2)), (Ala) 
x ch (l,2) = (T T n(l)n(2)) - (n(l))(n(2)), (Alb) 

where an abbreviation 1 = (ri,n) denotes a position ri 
and an imaginary time ti, and T T is the time ordering 
operator. Time dependence of a generic operator Q(r) 
is defined as Q(t,t) = e- TH Q(r)e rH . x sp(ch) (<7) in the 
main text is the Fourier transform of the above real-space 
representation. The two sum rules (Eqs. (2)), which play 
a central role in TPSC, originate from the definition at 
equal time, that is, X sp(ch) (l, 1+ ) with 1+ = ( r i. r i + s ) 



(S > 0). 

In the multi-orbital case, we consider correlation func- 
tions for S^l) = (c^(l), c^(l)) a (c Mt (l), Cl4 (l)) T , 

and n^(l) = 0^(1)0^(1) + 0^(1)0,4(1). These corre- 
lation functions for spin and charge channels are defined 
as 

X5EU(1,2) = <T T S* AM (1)SV(2)>, (A2a) 

X&(M) = {^5+^(1)5-^(2)), (A2b) 

Xtv € (1.2) = (T T n v (l)^(2)) - (n A/i (l))(n c „(2)). 

(A2c) 

The sum rules of Eqs. (8) and (9) come from the fol- 
lowing definitions at equal time, 



x5EU(M + ) 



xft*(M + ) 



((4 t (i)c M (i) - 4 4 (i)c^(i))(4 t (i+) C , t (i+) - 4 ; (i+)c H (i+))) 

(( C t t (i) C/lt (i) - ci i (i)c M (i)))(( C ^(i + ) C , t (i+) - c | i (i+) Ci4 (i + ))) 

2(o+ f (l)o /4 (l)o^(l+)o i/t (l+)) 

((4 t (l)c Mt (l) + ci i (l)c M )(4 r (l+)c yt (l+) + 4(l+)c^(l+))) 

((4 t (i) C/lt (i) + c t 4 (i) CM )(4 T (i + )c„ t (i + ) + c | i (i + ) c ^(i+))}. 



(A3a) 
(A3b) 
(A3c) 



r 



Appendix B: Ansatz for effective interactions in 
single-orbital case 

Following Vilk and Tremblay, 17 let us derive the 
ansatz, Eq. (4), used in TPSC calculations. The four- 
point vertex function, T aa -' , between electrons with spin 
a and a' is given by 



,5(1 -3)6(2 -4)5(2-1+) = 



3S ff (l,2) 
6G a ,(3,4Y 



(Bl) 



where G (T (1,2) and E CT (1,2) is the dressed Green's func- 
tion and self energy with spin a. The equation of motion 
and the Dyson equation leads to the relation, 



£ CT (1,1)G CT (1,2) 

= -U(T r [ct (l++)c ff (l + K(l)4 (2)]), 



with a = —a. Here a bar over a number means the in- 
tegral over position and imaginary time. The four-point 
correlation function can be approximated by the local 
correlation function and the Green's function as follows, 



C/(T T [4(1 ++ ) C ,(1 + ) C .(1)4(2)]> 

(%(iK(i)} 



u 



G ff (l,l+)G CT (1,2). 



(B3) 



(B2) By substituting Eq. (B3) into Eq. (Bl), we can obtain 



r ffff ,tf(i-3)a(2-4)<s(2-r 



^(1,2) f [tf^GWl, 1+)«(1- 2) 



<5G^(3,4) 
5 



<"t><™4-> 



<5G CT -(3,4) 



G,(l, !+)*(! - 2) + U , { W\f' { \'K} 6(l 2). 



(n t )( ni ) JG ff /(3,4) 



(B4) 



r 



The last term of this equation is proportional to b S(J i via and then contributes to the spin channel, U sp = T aS — 
SG-Cl 1 + ) ^aa- This just provides Eq. (4) for U sp - 



<JG ff /(3,4) 



<W<K1-3)J(4-1+), 



(B5) 
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Appendix C: Ansatz in multi-orbital systems point vertex function has orbital indices, A, fx, £ besides 

spin index, a. We here consider r^)^), which can 
In this section, let us extend the above-mentioned be written by only orbital-diagonal components, 
ansatz into the multi-orbital case. In this case, the four- 



r Mi ^)6(l - 3)6(2 - 4)5(2 - 1+) = 



6 



jj { n ficr n ficr) q 



as MMg (l,2) _ 5[S^(l,5)[G(5,6)G- 1 (6,2)] jxjxg ] 
6G (3,4) 6G vva >{3,4) 



^G ffff (l, 1+) - £(t/' - J) . (n ^ g) . G« g (l, 1+) ] 

s) (v)^) y 



-J7 



{n^a){nis) 8G vva >(Z,4) 



(n^in^) 8G VUC r'{3A) 



(CI) 



r 



Here, following the single-orbital case, we have intro- 
duced the following approximations, 



^flficrG fj.fj.crG ~ £/ 



{^fjcr^ficr) 

(nfjcr){nfjcf) 



Gfj,fj,cf(l, 1 + ), 



^fjfja^-r^cr^r^cr ~ 7 w ~ T ^ ^ ' ' ' 



1 V'ficr 



L fjcr 



G^(l,l+). 



Eq. (10a) for U?v„„ can be obtained from the first term 
of Eq. (CI), since the intra-orbital Coulomb interaction 
for the orbital fi is proportional to i^^j. In the same 
way, Eq. (10b) for [7^ and Eq. (10c) for U^ v - J s p 
with \x ^ v can be obtained from the second and the 
third terms in Eq. (CI), respectively. 
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